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Abstract. The orthogonal matching pursuit (OMP) is a greedy algorithm to 
solve sparse approximation problems. Sufficient conditions for exact recovery are 
known with and without noise. In this paper we investigate the applicability of 
the OMP for the solution of ill-posed inverse problems in general and in particular 
for two deconvolution examples from mass spectrometry and digital holography 
respectively. 

In sparse approximation problems one often has to deal with the problem of 
redundancy of a dictionary, i.e. the atoms are not linearly independent. However, 
one expects them to be approximatively orthogonal and this is quantified by the 
so-called incoherence. This idea cannot be transfered to ill-posed inverse problems 
since here the atoms are typically far from orthogonal: The ill-jjosedness of the 
operator causes that the correlation of two distinct atoms probably gets huge, 
i.e. that two atoms can look much alike. Therefore one needs conditions which 
take the structure of the problem into account and work without the concept of 
coherence. In this paper we develop results for exact recovery of the support of 
noisy signals. In the two examples in mass spectrometry and digital holography 
we show that our results lead to practically relevant estimates such that one may 
check a priori if the experimental setup guarantees exact deconvolution with OMP. 
Especially in the example from digital holography our analysis may be regarded 
as a first step to calculate the resolution power of droplet holography. 



AMS classification scheme numbers: 65J20, 94A12, 47A52 
1. Introduction 

We consider linear inverse problems, i.e. we are given a bounded, injective, linear 
operator K : B ^ H mapping from a Banach space B into a Hilbert space H. 
Moreover, we assume that for an unknown v £igK we are given a noisy observation 
with ||ti — ii^ll < e and try to reconstruct the solution of 

Ku^v (1) 
* Author to whom correspondence shall be addressed. 
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from the knowledge of v^. We are particularly interested in the case where the 
unknown solution u may be expressed sparsely in a known dictionary i.e. wc consider 
that there is a family £ := {ejjiez C B of unit-normed vectors which span the space 
in which we expect the solution and which we call dictionary. With sparse we mean 
here that there exists a finite decomposition of u with N atoms Si G E, 

u = ^^aiei with a, e R, ||q:||^o =: A'' < oo. 

In the following we denote with / the support of a, i.e. / = {i € Z | a, ^ 0}. For any 
subset J we denote £( J) := {e, | i S J}. 

This setting appears in several signal processing problems, e.g. in mass 
spectrometry [22] where the signal is modeled as a sum of Dirac peaks (so-called 
impulse trains): 

u = ttj 6{- — Xi). 

Other applications for instance can be found in astronomical signal prociessing 
problems or digital holography, cf. [40], where images arise as superposition of 
characteristic functions of balls with different centers x, and radii rj, 

ijez 

Typically K docs not have a continuous inverse and licncc. the solution of the operator 
equation (1) docs not depend continuously on the data. This turns out to be a 
challenge for the case where only noisy data with noise level ||w — < e arc 
available — as it is always the case in praxis. First a small perturbation e can cause an 
arbitrarily large error in the reconstruction u of "Ku = w^" and second no solution u 
exists if w*^ is not in the range of K. 

Inverse problems formulated in Banach spaces have been of recent interest and 
there axe a several results which deal with solving inverse problems formulated 
in Banach spaces, e.g. results concerning error estimates [1,8,17,20,27,28,36] or 
Landweber-like iterations or minimization methods for Tikhonov functionals, see 
e.g. [2-5,9,16,37,38]. 

In the following, an approximate solution of "Ku = v'^^' shall be found by deriving 
iteratively the correlation between the residual and the unit-normed atoms of the 
dictionary 

D:={di}ie^:={j^} . 

Note that since the operator K is injective we get that Kci ^ for alii gZ and hence 
the dictionary 'D is well defined. In any step wc select that unit-normed atom from the 
dictionary D which is mostly correlated with the residual, hence the name "greedy" 
method. To stabilize the solution of "Ku = v^" the iteration has to be stopped early 
enough. 

For solving the operator equation (1) with noiseless data and the case where 
only noisy data with noise-bound \\v — v^\\ < s are available we use the orthogonal 

matching pursuit, first proposed in the signal processing context by Davis ct al. in [30] 
and Pati et al. in [35] as an improvement upon the matching pursuit algorithm [31]: 
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Algorithm 1.1 Orthogonal Matching Pursuit 

Set fc and /° := 0. Initialize r" := (resp. r^ := v ioi e = 0) and vP := 0. 
while ||r''|| > e (resp. \\r''\\ ^ 0) do 
k:=k + l, 

ik e argsup {|(r''-i,di)| | di e D}, 
jfc :=/fe-iu{i4, 

Project u onto span£(J*^) 

u*^ := argmin {\\v^ — KuW^ \ u e span £(7'^)}, 
j,fe ._ _ Xu''. 

end while 



Remark that in infinite dimensional Hilbcrt spaces the supremum 

sup{|(r'="\d,)||d, eD} (2) 

docs not have to be realized. Because of that OMP has a variant — called weak 
orthogonal matching pursuit (WOMP) which does not choose the optimal atom in 
the sense of (2) but only one that is nearly optimal, i.e. for some fixed oj G (0, 1] it 
chooses some ik & ^ with 

|(r'=-\d,,)| >a;sup{|(r'=-i,di)|M,; e ©}. 

In [42] a sufficient condition for exact recovery with algorithm 1.1 is derived, and 
in [10] it is transfered to noisy signals with the concept of coherence, which quantifies 
the magnitude of redundancy. This idea cannot be transferee! to ill-posed inverse 
problems directly since the operator typically causes that the correlation of two distinct 
atoms probably gets huge. Therefore in [11,15] the authors derive a recovery condition 
which works without the concept of coherence. For a comprehensive presentation of 
OMP cf. e.g. [29]. 

The paper is organized as follows. In section 2 we refiect the conditions for exact 
recovery for OMP derived in [42] and [11, 15] and rewrite them in the context of 
infinite-dimensional inverse problems. Section 3 contains the main theoretical results 
of the paper, namely the generalization of these results to noisy signals. In section 4 
we apply the deduced recovery conditions in the presence of noise to an example 
from mass spectrometry. Here, the data are given as sums of Dirac peaks convolved 
with a Gaussian kernel. To the end of this section we utilize the deduced condition 
for simulated data of an isotope pattern. Another example from digital holography 
is concerned in section 5. The data are given as sums of characteristic functions 
convolved with a Presnel function. This turns out to be a challenge because the 
convolution kernel oscillates. Similar to section 4 we apply the theoretical condition 
to simulated data, namely to digital holograms of particles. The two examples from 
mass spectrometry and digital holography illustrate that our conditions for exact 
recovery lead to prac;tic;ally relevant estimates suc;li that one may check a priori if 
the experimental setup guarantees exact deconvolution with OMP. Especially in the 
example from digital holography our analysis may be regarded as a first step to 
calculate the resolution power of droplet holography. 

2. Exact Recovery Conditions 

In [42], Tropp gives a suSicient and necessary condition for exact recovery with OMP. 
Next, we list this result in the language of infinite-dimensional inverse problems. 
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Define the linear continuous synthesis operator for the dictionary D = {di} = 
{Kei/\\Kei\\} via 

D: e ^H, 

Since D is linear and bounded, the Banach space adjoint operator 

D* : H {fy = 

exists and arises as 

D*v = i{v,d,)),^^ = {{v,j^^))^^^. 

Note that the use of £^ and its dual arises naturally in this context. Furthermore, 
for J C Z we denote with Pj : — > £^ the projection onto J and with the 
pseudoinverse operator of A. With this notation we state the following theorem. 

Theorem 2.1 (Tropp [42]). Let a G £^ with suppa = I, u — source 
and V = Ku the measured signal. If the operator K . B ^ H and the dictionary 
£ = {ei}i^z fulfill the Exact Recovery Condition (ERC) 

sup ||(£>P,)td||,i <1, (3) 

then OMP with its parameter e set to recovers a exactly. 

Theorem 2.1 gives a sufficient condition for exact recovery with OMP. In [42] 
Tropp shows that condition (3) is even necessary, in the sense that if 

sup ||(DP7)^rf||^i > 1, 

then there exists a signal with support / for which OMP does not recover a with 
V = Ku = OLiCi. 

The ERC (3) is hard to evaluate. Therefore Dossal and Mallat [11] and Gribonval 
and Nielsen [15] derive a weaker sufficient but not necessary recovery condition that 
depends on inner products of the dictionary atoms of D{I) and T){fi) only. 

Proposition 2.2 (Dossal and Mallat [11], Gribonval and Nielsen [15]). Let a G £° 
with suppa = I, u — '^i^i^OiiCi he the source and v = Ku the measured signal. If the 
operator K : B ^ H and the dictionary £ = {ei}i^z fulfill the Neumann ERC 

sup^](di,dj)| + sup^|(di,rfj)] < 1, (4) 
then OMP with its parameter e set to recovers a. 



The proof uses a Neumann series estimate for PiD*DPi — this clarifies the term 
"Neumann" ERC. The proof is contained in [15]. 
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Remark 2.3. Obviously the Neumann ERG (4) is not necessary for exact recovery. 

A demonstrative example can be found in witli tlie signal v — (1,1,1,0)^ and 
the unit-normed dictionary D = {di := (1,0,0,0)^,^2 := 2-i/2(i^ q, 0)^, da := 
2-1/2(1, o,l,0)"^,rf4 := (0,0,0,1)"^}. Here with 7 = {1, 2, 3} and /"^ = {4} we get 

|(rfl,d4)| + 1(^2,^4)1 + 1(4,^4)1 = 

but 

\{d,,d2)\ + \{d,,d3)\ = V2>l, 

hence the Neumann ERG is not fulfilled. The ERG (3) is nevertheless fulfilled since 

in that case 1| (_DP/)'''(i4||£i = 0. OMP will then recover exactly, as one could expect 
by considering that just {di, (^2, 0(3} span the M^. 

This counter-example may be generalized by considering / c Z such that 

sup = and sup^^^ \ {di,dj)\ > 1. 

Here the Neumann ERG fails but for any signal with support / OMP will recover 
exactly since the atoms di, i e I, and dj, j G are uncorrelated and OMP never 
chooses an atom twice. 



Remark 2.4. The sufficient conditions for WOMP with weakness parameter u) e 
(0, 1] are 

sup \\{DPi)U\Ui <io 
deD(/(!) 

and ^ 

sup V dj)\ + - sup V I {di, dj)\ < 1, 

according to theorem 2.1 and proposition 2.2, respectively. They are proved 
analogously to the OMP case — same as all other following WOMP results. 

Usually for sparse approximation problems the behavior of the dictionary is 

characterized as follows. 



Definition 2.5. Let 9" := {,fi}iez be a dictionary. Then the corresponding coherence 
parameter /x and cumulative coherence fii{m) for a positive integer m are defined as 

M := swp\{fi,fj}\ 

and 

lii{m):= sup sup 

ACZ i^A . . 
|A|=m 

respectively. Note that /Ui(l) = /u and /xi(rn) < mfj, for all m € N. 

Since sup^^j X^jg/j^j IR. < A*i(^ - 1) and sup^^^c Eje/ IK.c^j)! < Mi(^) 
we get another condition in terms of the cumulative coherence, which is even weaker 
than the Neumann ERG: 
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Proposition 2.6 (Tropp [42]). Let a G with suppa = I, u = J2iei,^i^i 
source and v = Ku the m,6asur6d signal. If the operator K . B ^ H and the dictionary 
£ = {ei}i^z lead to a dictionary D which fulfills the condition 

Mi(Af-l) + /xi(Ar) < 1, (5) 

then OMP with its parameter s set to recovers a. 

Remark, that the condition in proposition 2.6 for ill-posed inverse problems 
might be unsuitable, since the typically compact operator causes that the coherence 
parameter fj, probably is close to one. Therefore the cumulative coherence can grow 
large with increasing support. 

Remark 2.7. Another major approach for solving sparse approximation problems is 
the basis pursuit (BP). Here one solves the convex optimization problem 

min \\a\\ii subject to -ftrY'ajei = v. 

This idea is closely related to Tikhonov regularization with sparsity constraint. Here 
the basic idea is to minimize least squares with £^-penalty, 

min \\KJ2oiiei - v\\% + 'y\\a\\ii. 

In [42] it is shown that the ERC (3) also ensures the exact recovery by means of BP. 
Since the proposition 2.2 and proposition 2.6 arc estimates for the ERC (3) and the 
proofs do not take into account any properties of the OMP algorithm these results 
hold here, too. 

3. Exact Recovery in the Presence of Noise 

In [10], Donoho, Elad and Temlyakov transfer Tropp's result [42] to noisy signals. 
They derive a condition for exact recovery in terms of the coherence parameter /x 
of a dictionary. This condition is — ^just as remarked in [10] — an obvious weaker 
condition. As already mentioned, in particular for ill-posed problems this condition is 
too restrictive. In the following we will give exact recovery conditions in the presence 
of noise which are closer to the results of theorem 2.1 and proposition 2.2. 
Assume that instead of exact data v = Ku G H only a noisy version 

= V + ri = Ku + T] 

with noise level \\v — v^\\ = \\r]\\ < e can be observed. Now, the OMP has to stop 
as soon as the representation error r'^ is smaller or equal to the noise level e, i.e. if 

Theorem 3.1 (ERC in the Presence of Noise). Let a G £° with suppa = /. Let 
u = ^i^i source and = Ku + rj the noisy data with noise level \\r]\\ < s 

and noise-to-signal-ratio 

snp \{r],di)\ 
iez 

' mm aj Aej 
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// the operator K : B ^ H and the dictionary £ = {ei}i£-L fulfill the Exact Recovery 
Condition in Presence of Noise (eERC) 

sup ||(/?P,)t4,x <l-2r,/„— — — , (6) 

iei jei 

and sup Yl lidijdj)] < 1, then OMP recovers the support I of a exactly, 
iei jei 

Proof. We prove the eERC analogously to [42, theorem 3.1] by induction. Assume 
that OMP recovered the correct patterns in the first k steps, i.e. 

iei'' 

with l'^ c /. Then we get for the residual 

r'^ := - Kv!' = v + r]- Ku'^ = K(^^{ai - a'l)e^ + 1] 

= ^\\Kei\\{ai-a^)di + T], 
iei 

hence the noiseless residual s*^ := J2 \\-^^i\\i^i ~ cei)di has support /. The correlation 
\{r'',di)\, i e Z, can be estimated from below and above respectively via 

|(r^d,)| = |(s'= + 77,d,)|||(s^d,)|T|(ry,di)|. 

Hence with 

sup \{r'',d^)\ < sup \ {s'',di)\+ sup \{r],di)\ 
iei° iei° »ez 

and 

sup I (s^ di) I - sup I (r?, d,) I < sup | (r^ di) | 
iei iei iei 

we get the condition 

||Pjcr'*s*^||^=o +sup|(j7,di)| < \PiD*s%io. -sup|(77,rfi)|, 

which ensures a right choice in the (fc + l)-th step. Since {PiD*y PjD* is the 
orthogonal projection onto and supps*' = / we can write 

5^= = {PjD*yPiD*s''. 

With this identity we get the sufficient condition for OMP in presence of noise 

\\P,oD*iPrD*VP,D*s%^ 



\\PiD*s>'\\eoo \\PiD*s>'\\e^ 
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Consequently, since |lP^ci:'*(P/i:'*)1'P/D*s'=||^=o < \\Pj:,D* {PiD*y\\i^ \\PiD*s%o. 
and with the definition of the adjoint operator D* , we get the equivalent sufficient con- 
ditions for a correct choice in the {k + l)-th step 



sup I (f/. (■/,) 



\\PiD*s^\\eo 
Obviously, on the one hand we get 

sup \\{DPi)Ui\\e. < sup ||pPz)t^A(i,||,i = \\{DPi)^DPj,\ 



supp 13=1 



and on the other hand, since {DPiY is linear, we get 

sup ||(DP,)t VArfijj^i < sup V|/?i| sup ||(£'P7)^rf|Ui 

supp/3=/ supp/3— / 

= sup \\{DPj)U\\e. 
deii(/(i) 

This shows that 

sup I (77,(1,) I 

||(£>Pj)t£)P^o 11^1,^1 = sup ||(DP7)td||^i < 1 - 2 ■ 



uiij^ V-^-'j; ^ ^ II n 7-1* HI 

de©(/C) ||P/D*s'=||^oc 

is another equivalent condition for exact recovery. The last thing we have to afford to 
finish the proof is an estimation for the term ||P/£)*s'^||foo from below. 

In the first step this is easy, since r^ = resp. s^ = v with v = Ku = 
KJ2iei cti^i- With that we get 

||P/£>*s°||^oo = ||P,D*t;||^oc =supKt;,d,)| =sup|^ai||irei||(di,(i 



> 



Y,(^i\\Kei\\{di,di) > |anil^ei||(l-sup^|(di,rf,)|) 

for all Z e /, hence in particular 

||P/I?*f||^oo >min|ai|||i^e,||(l-supV|(d„d,)|). 

To prove this for general k we successively apply this estimation. Here, again, we 



get 



\\PiD*s''\\^ = sup\{s\dj)\ = sup\j2{a'^ - ai)\\Keim,dj)\ 



jei jei 



> 



Y^i^i - ai)\\Kei\\{di,di)\ > \ai\ ||i^en|(l - sup^ Krfi,rf,)|) 



iei jei 
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for all Z e /, Z ^ Z*', hence, since OMP never chooses an atom twice, in particular 

\\PiD*s''\\eo. > min|ai| ||Kei||(l - sup V 



□ 



In particular, to ensure the eERC (6) one has necessarily for the noise-to-signal- 
ratio r^/a < 1/2. For a small noise-to-signal-ratio the eERC (6) approximates the 
ERG (3). A rough upper bound for supj^^ liv^dt)] is e and hence, one may use 



re/a 



< 



min lajl llifej 



Similar to the noiseless case, the eERC (6) is hard to evaluate. Analogously to 
section 2 we now give a weaker sufficient recovery condition that depends on inner 
products of the dictionary atoms. It is proved analogously to proposition 2.2. 

Proposition 3.2 (Neumann ERG in the Presence of Noise). Let a € with 
suppa = /. Let u = J2iez ^i^i source and v'^ = Ku + ij the noisy data with 

noise level \\r]\\ < e and noise-to-signal-ratio r^/^ < 1/2. // the operator K : B ^ H 
and the dictionary £ = {ei}i^z fulfill the Neumann eERG 



sup ^\{di,dj)\-\- sup ^\{di,dj)\ <l-2 r^/a, (7) 
then OMP recovers the support I of a exactly. 



Remark 3.3. The according suffient conditions for exact recovery with WOMP with 
weakness parameter w e (0,1] for the case of noisy data with noise-to-signal-ratio 
re/a < w/2 are 

^^SUp \iiDP,)U\U. < . -^re/a^^^^^^^J^^, 



and 



2r 



e/a 



sup V|(di,dj)| + - sup V|(di,dj)| < 1 

analog to theorem 3.1 and proposition 3.2, respectively. 

Same as for the noiseless case we can give another even weaker condition in terms 
of the cumulative coherence: 

Proposition 3.4. Let a G with suppa — I. Let u = '^i^iOiiCi he the source 
and v'^ = Ku -\- rj the noisy data with noise level ||ry|| < e and noise-to-signal-ratio 
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fs/a < 1/2- If the operator K : B ^ H and the dictionary £ = {ei}i^z lead to 
dictionary D which fulfills the condition 

/zi(Ar-l) + /zi(Ar) <l-2r,/„, (8) 

then OMP recovers the support I of a exactly. 

Remark that theorem 3.1 and proposition 3.2 just ensure the correct support I. 
The following simple proposition shows that the reconstruction error is of the order 
of the noise level. 

Proposition 3.5 (Error bounds for OMP in presence of noise). Let a G £^ with 

suppa ~ I. Let u = X^iez'^^'^* source and the noisy data with noise level 

\\v'^ — t;|| < £ and noise-to-signal-ratio r^/cx < 1/2- If the eERC is fulfilled then there 
exists a constant C > such that for the approximative solution a. determined by OMP 
the following error bound holds, 

\\a — a\\(i < Ce. 

Proof. Since the eERC is fulfilled OMP recovered the correct support /, i.e. 

a = argmin{||t;^ - J2iei \ a€f{I)}. 

With the help of the operator A : i^{I) H defined by Aa = J2iei '^i^^i 
equivalently written as 

A*Aa = A*v^. 
Note that A* A : £'^{I) f{I) is just the matrix 

{A*A)ij = {Kei,Kej). 

For the error we get 

lia - a||,i = \\A^v' - v)\\er < \\A^\\H,e^ \\{v' - v)\\h = Ce. 

□ 

Remark 3.6. We remark again on an exact recovery condition for BP. Unlike the 
section 2 where the results can be transfered to BP, see remark 2.7, this is not possible 
for the presence of noise: To prove theorem 3.1 we used properties of the OMP 
algorithm which are not valid for BP. 

For the case of noisy data in [10] an exact recovery condition for BP is derived. 
This condition depends on the coherence parameter /i. Since in this paper the focus 
is on the greedy solution of inverse problems we give up on deriving stronger results 
for BP which are closer to the results of this section. 

4. Resolution Bounds for Mass Spectrometry 

Granted, to apply the Neumann conditions of proposition 2.2 and proposition 3.2, 
respectively, one has to know the support /. In this case there would be no need to 
apply OMP — one may just solve the restricted least squares problem, i.e. project onto 
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span £ (7). For deconvolution problems, however, with certain prior knowledge the 
Neumann ERC (4) rcsp. Neumann eERC (7) arc easier to evaluate than the ERC (3) 
resp. eERC (6) especially when the support / is not known exactly. In the following 
we will use the weaker conditions exemplarily with impulse trains convolved with 
Gaussian kernel as e.g. occurs in mass spectrometry, cf. [22]. 

4-1- Analysis 

In mass spectrometry the source u is given — after simplification — as sum of Dirac 
peaks at integer positions i gZ, 



with |suppQ;| = |/| = A''. Since the measuring procedure is influenced by Gaussian 
noise the measured data can be modeled by a convolution operator K with Gaussian 
kernel ^ 

i.e. the operator under consideration is Ku = k * u. As Banach space B we 
may use the space A4 of regular Borel measures on M (which contains impulse 

trains if the coefficients are summable) and as Hilbert space H the space L^(IR). 
We form the dictionary £ of Dirac peaics at integer positions and hence, we have 
2) = {6{- — i) * k} = {k{- — i)}, since \\k{- — i)\\L2 = 1, i gZ. 

To verify the Neumann ERC (4) and Neumann eERC (7) respectively, we need 
the autocorrelation of two atoms k{- — i) and k{- — j). In L^(M,M) it arises as 

{n{-i),n{-j))L.= j -4-exp(-(^)exp(-(^) da. = exp (- (^), (10) 

R 

which is positive and monotonically decreasing in the distance |i— i|. If we additionally 
assume that the peaks of any source u have the minimal distance 

p := min |i — j|, 

i,iGsupp a 

then w.l.o.g. we can estimate the sums of correlations from above as follows. For 
p e N we get for the correlations of support atoms 

1JV/2J VN/2\ 

sup^|(d„d,)|<2 ^ («,«(. -jp)) = 2 ^ exp(-(|f)^). 
jei j=i j=i 

For the correlations of support atoms and non-support atoms we have to distinguish 
between two cases for p. For p > 2 we get 

sup^\{di,dj)\ < sup {k{- -i),K{- - jp)) = sup ^ exp ( - 
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and for p = 1 

LAf/2j+l [N/2i+l 

sup^Krf,,rf,)| <2 J2 («,«(--i))=2 Yl exp(-3j). 

With that we can formulate the Neumann ERC and the Neumann eERC for Dirac 
peaks convolved with Gaussian kernel. 

Proposition 4.1. An estimation from above for the ERC (i.e. r^/^ o,nd eERC 
(i.e. < r^/a < ^) for Dirac peaks convolved with Gaussian kernel is for p>2 

VN/2\ . 2 VN/2\ 2 

i.JP) \ , „ V- / {i-Jp) ' 



and for p = 1 

2Ee-p(-^)+2 E ^-p(-£2)<l-2re/a. 

This means that we are able to recover the support of the impulse train with OMP 
exactly from the convolved data if the above conditions are fulfilled. 



lN/2\ ^ .2 , LJV/2J+1 



2 LJ'/^J^i -2 



Remark 4.2. Remark that the case p = 1 of proposition 4.1 coincides more or 
less with the recovery condition in terms of the cumulative coherence, since for odd 
N we get Hi{N) = 2 ^fIlexp{—j-^/{Aa^)). Summing up just over a subset of 
pZ := {j G Z\j/p G is not a feasible estimation, since for the support I we 
allow any point i G Z and not only atoms of the sub-dictionary ©(pZ). This turns 
out to be the main disadvantage of the coherence condition: It does not distinguish 
between support and non-support atoms and is hence in some applications a clearly 
weaker estimation. 



Remark 4.3. If the cardinality of the support N is unknown one could replace the 
finite sums by infinite sums. Obviously these sums exist since the geometric scries is 
a majorizing series. With l representing the imaginary unit they can be expressed in 
terms of the Jacobi theta function of the third kind, l^3{z, q) := ^^^=-00 1'' exp(2ji,2;). 

The condition of proposition 4.1 is plotted for some combinations of a, p and rg/Q, 
with unknown N in figure 1. The colored areas describe the combinations where the 
Neumann ERC is fulfilled. 

Often for dcconvolution problems the autocorrelation of two atoms \ {d{- — i), d{- — 
is not monotonically decreasing in the distance and it obviously depends on 

the kernel k. However, if the correlation of two atoms can be estimated from above 
via a monotonically decreasing function w.r.t. an appropriate distance then we can 
use a similar estimate. We do this exemplarily for an oscillating kernel in section 5 
namely, for Fresnel-convolved characteristic functions as appear in digital holography. 
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Figure 1. eERC for combinations of cr, p and r^i^. 

Remark 4.4. Wc remark on a possible fully continuous formulation of OMR We 
assume that we are given some data 

V = U = Ui k{- — Xi) 

and that we do not know the positions Xj. We allow our dictionary to be uncountable, 
i.e. we search for peaks at every real number. Note that here i G Z does not represent 
the set of possible positions for peaks but it is an index set for continuous positions 
Xi G M. 

In the first step of the matching pursuit we correlate v with k{- — x) and take 

that X which gives inaxiinal c;orrelatioii. In the spec;ial c;ase of the Gausssian blurring 
kernel (9) this amounts in finding the maximum of the function 

f{x) = \{v, k{- - a;))| = ^ ai{K{- - Xi),K{- - x)). 

From (10) we see that this is 

/(x) = ^ exp ( - . 
iez 

It is clear that any maximum of / is unlikely to be precisely at some of the XiS, albeit 
very close. A detailed study of this effect goes beyond the scope of this paper and we 
present a simple example. 

Let us assume that we have two peaks, one at and one at x\: 



u = aoS{-) + ai6{- — xi). 



(11) 




Figure 2. Error of the first step of the matching pursuit for the signal (11) with 
ao = 2, ai = 1 and xi = 1. The variable a is the variance of the Gaussian kernel 
and Q is the position at which the matching pursuit locates the first peak. 



Moreover, we assume that ao > ai, i.e. the peak in zero is higher. The matching 
pursuit will find the first peak at the maximum of the function 

(x'^ \ ( ix — X\)^ \ 

and hence at some root of 

= -2^2 (aoa:exp(- ^) + ai(a; - a;i) exp ( - ^"^4^2^ ))• 

The error, that the matching pursuit makes is hence the error q in the root of /' near 
zero. In figure 2 it is shown, how the root of /' close to zero depends on the variance 
a. One observes that the error q is smaller than the variance a by some orders of 
magnitude. 

As a final remark we mention that we measured the error not in some norm but 
only the distance of the 5-peaks. This corresponds to the so called Prokhorov-metric 
which is a metrization for the weak-* convergence in measure space. 



Numerical Examples 

We apply the Neumann eERC of proposition 4.1 to simulated data of an isotope 
pattern. Here the data consist of equidistant peaks with different heights. In our 
example we use four peaks with a distance of p = 5 and heights of 130, 220, 180 
and 90, cf. the balls at the top of figure 3. After convolving with Gaussian kernel 
with a = 1.125 we apply a Poisson noise model. This is realistic, because in mass 
spectrometry a finite number of particles is counted. 

In the first example with low noise (mean and variance of 1.5 for regions without 
peaks) the Neumann eERC is fulfilled and hence OMP recovered the support exactly, 
see middle of figure 3. However, the condition is restrictive: For the second example 
the signal is disturbed with huge noise (mean and variance of 30 for regions without 
peaks) and the Neumann eERC is not fulfilled. Certainly, OMP recovered the support 
exactly, see bottom of figure 3. 
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Figure 3. Simulated isotope pattern. Top; support and Gaussian-convolved 
data without noise. Middle; low noise, Neumann eERC satisfied. Bottom; high 
noise, Neumann eERC not satisfied but still exact recovery possible. 
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5. Resolution Bounds for Digital Holography 

In digital holography, the data correspond to the diffraction patterns of the objects [12, 
24]. Under Presnel's approximation, diffraction can be modeled by a convolution with 
a "chirp" kernel. In the context of holograms of particles [18,19,44], the objects can be 
considered opaque (i.e., binary) and the hologram recorded on the camera corresponds 
to the convolution of disks with Fresnel's chirp kernels. The measurement of particle 
size and location therefore amounts to an inverse problem [39,40]. 



Analysis 

We consider the case of spherical particles which is of significant interest in applications 
such as fluid mechanics [32,45]. Wc model the particles j E {1, . . . , N} as opaque disks 
Br{- — Xj,- — yj,- — Zj) with center {xj,yj,Zj) € M"^, radius r and disk orientation 
orthogonal to the optical axis {Oz). Hence the source u is given as a sum of 
characteristic functions 



N 



N 



Xj , • yj,- 



i) =■ 



The real values are amplitude factors of the diffraction pattern that in praxis depend 
on experimental parameters, cf. [40,43]. 

To an incident laser beam of (complex) amplitude Aq and wavelength A the 
amplitude A in the observation plane, i.e. at depth z = 0, is well modeled by a 
bidimensional convolution ® w.r.t. {x,y). In the following l represents the imaginary 
unit. Then, with 5xj,yj denoting Dirac's peak located at {xj,yj) and hzj the Presnel 
function, 

hzj {x, y) = exp (^i^E?^ , with E? := x^ + y'^, 
the amplitude A : — > C arises as 



N 



A = Aq 



{Xj®hz^®5x.,y^) 



Remark that hz^ ® 5^^ denotes the shifted Presnel function. 

One difficulty occurring at digital holography inverse problems is that in praxis 
only the absolute value of A can be measured by the detector and the phase gets lost. 
The measured intensity consequently arises as 



N 



G = \Af = l^oj' [l - 2 ^ aj {xj ® Mhzj) ® 

TV TV 

+ ^^oti{xi® ® 5x,,yi)aj{xj ® h_zj ® Sxj,yj) 



--1 ]=i 



Since the second term is dominant over the third one for x small, the intensity is 
classically linearized [40,43]: 



JV 



l-2^aj {xj ® Be{hzj ) ® S 



(12) 
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Analogously to section 4 we will next derive the Neumann ERC and the Neumann 
ERC in presence of noise for the operator equation (12), ^ aj Xj '-^ G. Here for fixed 
{xj,yj,Zj) the associated (not necessarily unit-normed) atoms d^^ € D have the form, 

d^^{- -Xj,- -yj) ■■= XbA- -Xj,- -yj)®'R'e{hzi)®Sxj,vr (13) 

As before the first step is to calculate the norm of an atom and the correlation of two 
distinct ones. Therefore we need some properties of the Fresnel function. 

Proposition 5.1. For the convolution of the Fresnel function we have the 
properties [26] 

hzi ® = hzi+z2 , for all zi,Z2e K, 

hz ® h-z = S, for all 2: e K. 

With that and hz — h-z we get for the real part of the Fresnel function 

Re{hzJ ® Re^hz^) = KRe^hz^+z^) + Re{hzi-Z2)), for all zi,Z2 € K, 
Re{hz) ® Re{hz) = ^{5 + Re{h2z)), for all zgR. 

Another important property is that the convolution of a function with the Fresnel 
function — the so-called Fresnel transform — can be related to a direct multiplication 
with its Fourier transform which is defined by: 

^f{C,v) = j f{x,y) ex.p{-2Tn{x^ + yr])) dx dy. 

Proposition 5.2. Let f G L^{M?) and hz be a Fresnel function. Then 

(/ ® hz) ii, v) = .^{^Az hz /} ^) hz{£,, v). 

Proof. Let / G L^(IR^), z G M and h^ the corresponding Fresnel function. Then 
rearranging yields to the statement, 

{f®hz){^,n) = j /(x,y)^exp(^((a;-0' + (2/-r?)^)) dx dy 
= ri^exp(^(^^ + r?^))|/(a;,y)exp(^(a=^ + y^))exp(-2^.(|| + f^)) Ax dy 

R2 

= T{i\zKf}{l^,^)hz{i,ri). 

□ 

Remark 5.3. In praxis / has a bounded and small support w.r.t. V\z. With 

(x^ + ?/^)max denoting the maximal spatial dimension of / resp. the maximal spatial 

/ 2_i_ 2 \ 

extend of the corresponding particle the so-called far-field condition ^ 
holds in the proof of proposition 5.2, cf. [43]. In [40] e.g., particles of radius at about 
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50iim are illuminated with a red laser beam (wavelength 630nm) and distance to 
camera of about 250mm. Thus the term {t? + y'^)maxl{^z) « 3 • 10~^ and hence 
g-jj-p ^££(EJJLi^ is approximatively 1. Under the far-field condition we can estimate 



(14) 



With that for the complex valued difi^r action, with := ^'^ -\- rf and J^, denoting the 
first kind Bessel function of order v we get 



since J^XeXp) = 27rr^ 
real valued intensity atom we get 



Jl(27rrp) 
2'jrrp 



holds (Airy's pattern, vide infra). With that for a 



XBr ® Ile(/i2) = Re(xi3^ ® h;,) « '-^^y^^Pj (x^'^^j' 

which corresponds to the model given by Tyler and Thompson in [43]. 

Back to the correlation and — as a special case — the norm of an atom: The 
correlation appears as the autoconvolution, namely 



d-Zii' -^ii' yi)idzj{' ■^ji' Vj) I 

= / d^.{x,y) d^.{x - {xj -Xi),y- {yj - yi)) dx dy 



= {dzi ® d:,.){xj - Xi,yj - yi). 



In the following we assume that all particles are located in a plane parallel to the 
detector, i.e. z := Zi is constant for all i. Then the autoconvolution of an atom 
appears as _ _ 

d^®d^ = xb^ ® XBr ® Re(/i^) ® Re{h^). 



With proposition 5.1 and the formula 

'2r2cos-i (|;) - ^^/4r^^-p^2 for 4r'^ > p^ 
0, else. 



C{p) = {xBr®XBr){p) 



we get 



dz ® d. 



C + C® Re(/i2z) 



(15) 



6 + Re(/i2z) = 
With remark 5.3 and since J^C is real valued we get 

C ® Re{h2z) = Re(C ® h2z) « Re(j^C(-/Az) /la^) = :FC{-/Xz) Re(ft,2z) 

= ^XbA-/^z) ^XbA'/^z) Rc{h2z)- 

In physics it is well known that the Fourier transform of a disc is the Bessel cardinal 
function, Jinc(a;) := Ji{x)/x, since it is the diffraction of a circular aperture at infinite 
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distance. Nevertheless, for the sake of completeness and mathematical beauty we will 

illustrate this computation: Since the Fourier transform of a radial function is the 
Hankel transform of order zero (also known as Bessel transform of order zero), cf. [41, 
Theorem IV.3.3, page 155], the Fourier transform of appears, for := + rf, 
as 

r 2'Krp 

^XbM = 27ry 5Jo(27rp5) d5= ^ j S Jo{S) dS. 



In order to solve this definite integral we use J SJo{S) dS = SJi{S), cf. [21, equation 
5.52 1.], and get 

Ji(27rrp) 



2'Krp 



hence the Fourier transform of the circle-circle intersection C appears as 

With that result we can easily calculate the norm of an atom d^: Since C(0) = irr^, 
J^C{0) = {J'XBrWf = (/XB. da;)^ = ttV^ and /i2z(0) = we obtain 

= |4®4|(0) ~ iTrr^. 
Hence for fixed z we can represent the associated unit-normed atoms G D, with 



2 .. 



X + y , via 



dz ■■= 



'27rr 



\\dz 



(16) 



In figure 4 the centered atom for a particle of SOpm radius is displayed which is 
illuminated with a red laser beam (wavelength 630nm) in a distance of 250mm to the 
camera. 

The autoconvolution for general p and hence the correlation of two atoms 
dz{- —Xi, ■ —yi) and dz{- — Xj,-—yj) with distance distance p = {{xj—Xi)'^ + {yj—yi)'^)^ 
in digital holography emerges as 



dzi- - Xi, - -yi),d;,{- - Xj,- - yj)) = \dz®dz\{p) 



1 



\dz®dz\{p) 



1 

7rr^ 
C{p) 



C{p) + TC 



|Re(/i2.(p)) 



7rr 



1 ^il^-Kr 
2 +7-^i(— 



V \z 



(17) 



where sine denotes the normalized sine cardinal and is defined via sinc(a;) := 
sin(7ra;)/7ra;. 

The correlation in digital holography (17) is not as easily valuable as in mass 
spectrometry, because it is not monotonically decreasing in the distance p due to the 
oscillating Bessel and sine functions. To come to an estimate from above which is 
monotonically decreasing we use bounds for the absolute value of the Bessel functions 
J\. In [25] Landau gives estimates for | Ji/(a;)| for a; > and > 0, namely 



\J,(x)\ < m.in{bLV-^/^,CLX-^/^}, 



(18) 
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''-=^Ssf(-'-*(r) + -'*(5-'))==»-"^«- 

cl := sup a; 5 Jo (a;) w 0.7857. 

a;>0 

In addition sine cardinal obviously is bounded from above via 1 and 1/x and hence 
we have 



which now is monotonically decreasing in p. Figure 5 illustrates the oscillating part 
of the correlation (17) and its corresponding upper bound for two particles of SOpm 
radius which are illuminated with a red laser beam (wavelength 630nm) in a distance 
of 250mm to the camera. 



correlation 




-4- 10-4 



-8- 10-4 



(pm) 



Figure 4. Unit-normed, cen- 
tered atom of particles of ra- 
dius 50ijm, illuminated with 
a red laser beam (wavelength 
630nm) and distance to camera 
of 250mm. 



Figure 5. The oscillating part 
of the correlation of two atoms 
with distance p and its cor- 
responding monotonically de- 
creasing estimate (same set- 
tings as in figure 4). 



Remeirk 5.4. In [23] Krasikov gives more precise estimation for J^{x), namely, for 
u > -1/2, c := {2u + l){2v + 3) and x > ^^ + ^2/3/2, 

2 4{4x^-{2u + l)i2u + 5)) 

7r((4a;2-^)3/2_^) • 

With that (asymptotically \dz ® dz\{p) ~ p~^) instead of Landau's rough bound (18) 
(asymptotically \dz ® d^Kp) ~ p~^) one can get a more precise recovery condition 
for digital holography. Since this technical computation is beyond the scope of this 
theoretical paper we postpone it here. 
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With this estimation we will come to a resolution bound for droplets jet 

reconstruction, as e.g. used in [40]. Here monodispcrsc droplets, i.e. they have the 
same size, shape and mass, were generated and emitted on a strait line parallel to the 
detector plane. This configuration eases the computation of the Neumann ERC and 
the Neumann ERC in presence of noise. Analogously to mass spectrometry we define 
that the particles appear at some selected points i G AZ := {i e Z | ^ € Z}, where 
the parameter A describe the dictionary refinement. If we additionally assume that 
the particles have the minimal distance p € AN, then the sum of inner products of 
support atoms D{I) and non-support atoms can be estimated from above. For 

p > A we get 

LiV/2J 

supY,\{di,dj)\< Yl ^C(ip) + imin{6i,ci(^)i(ip)-i}min{l,^(ip)-2} 
jei 3=1 



and 



LJV/2J 

sup < sup J2 ^c-dip-*!) 

+ imin{6i,ci(^)i|ip-i|-i}min{l,^|ip-i|-2}. 

Proposition 5.5. An estimation from above for the ERC (i.e. r^/a = 0) and eERC 
(i.e. < r^ia < \) for characteristic functions convolved with the real part of the 
Fresnel kernel is for p > A 



LJV/2J 

E ^C{jp) + \ min {6i, ci(^) § (jp)"! } min {l, ^(jp)"^} 

+ si^P E ^C'(lip-iAi) 

l<i<2: 1^ 3 = -\NI1\ 

+ imin{6i,ci(^)i|ip-iA|-i}min{l,?A^|jp-iA|-2}| 



Remark 5.6. Same as before for mass spectrometry: If the cardinality of the support 
N is unknown one could replace the finite sums by infinite sums. These sums exist and 
can be expressed in terms of the Hurwitz zeta function C^(y, q) := X^°lo(9 + i)"' for 
f > 1, q > 0, and the Riemann zeta function ({v) := (^{u, 1) = X^^lii""^, respectively. 

The condition in proposition 5.5 seems not to be easy to handle. However, in 
praxis all parameters are known and one can compute a bound via approaching from 
large p. As soon as the sum is smaller than 1, it is guaranteed that OMP can recover 
exactly. A typical setting for digital holography of particles is the usage of a red 
laser of wavelength A = 0.6328pm and a distance of z = 200mm from the camera, 
cf. [40]. In figure 6 the condition of proposition 5.5 is plotted for particles with typical 
radii r e {5, 15, 25, 35, 50, 75}iim. In the computation the asymptotic formula is used. 
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i.e. for an unknown support cardinality N. For the dictionaries a corresponding 

refinement of A = r/2 was chosen. The colored areas describe the combinations 
where the Neumann ERC is fulfilled and hence OMP recovers exactly. 




400 600 800 1000 1200 1400 



Figure 6. sERC for combinations of p and r^/^ with corresponding dictionary 
refinement of A = r/2. For particles the radii r £ {5, 15, 25, 35, 50, 75}iim and 
the asymptotic formula (19) for an unknown support cardinality N are used. 

Numerical Examples 

We apply the Neumann eERC (7) to simulated data of droplets jets. For the simulation 
we use the same setting as above, i.e. a red laser of wavelength A = 0.6328iim and 
a distance of z = 200mm from the camera. The particles have a diameter of 100 
microns and the corresponding dictionary the refinement of 25iim. Those parameters 
correspond to that of the experimental setup used in [39, 40] 

After applying the digital holography model (12) we add Gaussian noise of 
different noise levels and in each case of zero mean. For the coefficients we choose 
2aj = 10 for all i G suppa. The figure 7 shows three simulated holograms with 
different distances p and noise-to-signal levels r^/Q,. For all three noisy examples 
in the right column all the particles were recovered exactly. However, only for the 
image on top{p ~ 721iim) the condition of proposition 5.5 holds. In the second 
image in the middle of the figure the particles have a too small distance to each 
other(/9 SSOpm) and even for the noiseless case the condition is not fulfilled. The 
last image (p w 721iim) was manipulated with unrealistically huge noise so that here 
the condition of proposition 5.5 is violated, too, cf. figure 6. 

6. Conclusion and Future Prospects 

In this paper we gave exact recovery conditions for the orthogonal matching pursuit for 
noisy signals that work without the concept of coherence. Our motivation was to treat 
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Figure 7. Simulated holograms of spherical particles. In the left column the 
noiseless signals v are displayed. For reconstruction the noisy signals of the 
right column are used. The dots correspond to the location of detected particles 
with OMP. The algortihm recovered all particles exactly, however, the condition 
of proposition 5.5 was just fulfilled for the image on top right. In the image in the 
middle the particles have a too small distance to each other and at the bottom 
the image was manipulated with unrealistically high noise. 
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ill-posed problems, and in particular, two problems of convolution type. We obtained 

results on exact recovery of the support for noiseless and even noisy data. Moreover, for 
noisy data there is a simple error bound in proposition 3.5 which shows a convergence 
rate of 0{e). The rate of convergence resembles what is known for sparsity-enforcing 
regularization with £p penalty term for < p < 1 [6,13,14], moreover, our results also 
guarantee the exact recovery of the support. 

In two real- world applications we showed that these condition lead to computable 
conditions and hence, are practically relevant. A main tool here was, that the atoms 
in the dictionary are shifted copies of the same shape and that the correlation of 
the atom depends on the distance of the atoms only. Once there is a sufficiently 
decaying upper bound for the correlation, we are able to apply the Neumann ERC (4) 
and the Neumann eERC (7) and obtain computable conditions for exact recovery as 
illustrated in the examples in section 4 and 5. However, experiments indicated that 
the conditions for exact recovery from theorem 2.2 and 3.2 are too restrictive. An idea 
to come to a tighter exact recovery condition is to bring in more prior knowledge, as 
e.g. a non-negativity constraint, cf. [7]. We postpone this idea for future work. For 
the particle digital holography application even more prior knowledge may be taken 
into account, since the objects are not only non-negative but also all apertures have 
the same denseness, i.e. aj is constant for all i G I. 

As discussed in remark 4.4, a straightforward generalization of our approach to 
fully continuous dictionaries runs into problems. Especially it seems that there is little 
hope to obtain exact recovery of the support, but maybe one may obtain bounds on 
how accurate the support is localized. This is strongly related to the structure of the 
dictionary (e.g. that is consists of shifts of the same object) and of course related to 
the correlations. 

Finally a further direction of research may be to investigate other types of pursuit 
algorithms like regularized orthogonal matching pursuit [34] or CoSAMP [33]. 
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